#include "sound-wave.hh"
#include <math.h>
namespace Dune {
  namespace Initial {
    SoundWave::SoundWave(Dune::ParameterTree& ptree){
      set_initcond_para("dim",1.0);
      set_initcond_para("nvar",1.0);
      set_initcond_para("rho_rectangle",0.0);
      set_initcond_para("rho_circle",0.0);
      set_initcond_para("rho_plane",0.0);    
      //
      set_initcond_para("wid1",0.0);
      set_initcond_para("wid2",0.0);
      set_initcond_para("wid3",0.0);
      set_initcond_para("r",0.0);
      set_initcond_para("k1",0.0);
      set_initcond_para("k2",0.0);
      set_initcond_para("k3",0.0);
      set_initcond_para("cen1",0.0);
      set_initcond_para("cen2",0.0);
      set_initcond_para("cen3",0.0);
      //
      set_initcond_paras(ptree);
    }

    SoundWave::~SoundWave(){
    }
    
    void SoundWave::set_initcond_paras(Dune::ParameterTree& ptree){
      para["rho_rectangle"]=ptree.get("initial.rho_rectangle",0.0);
      para["rho_circle"]=ptree.get("initial.rho_circle",0.0);
      para["rho_plane"]=ptree.get("initial.rho_plane",0.0);
      para["dim"]=ptree.get("grid.dimension",0.0);
      para["cen1"]=ptree.get("initial.cen1",0.0);
      para["cen2"]=ptree.get("initial.cen2",0.0);
      para["cen3"]=ptree.get("initial.cen3",0.0);
      para["r"]=ptree.get("initial.r",0.0);
      para["wid1"]=ptree.get("initial.wid1",0.0);
      para["wid2"]=ptree.get("initial.wid2",0.0);
      para["wid3"]=ptree.get("initial.wid3",0.0);
      para["k1"]=ptree.get("initial.k1",0.0);
      para["k2"]=ptree.get("initial.k2",0.0);
      para["k3"]=ptree.get("initial.k3",0.0);
      // DEBUG
      std::cout << "initial init: rho_rectangle="  
        << para["rho_rectangle"] << std::endl;
    }

    void SoundWave::initial_condition(double* location, double* holder){
      if (para["rho_rectangle"]==1.0)
        rho_rectangle(location,holder);
      if (para["rho_circle"]==1.0)
        rho_circle(location,holder);
      if (para["rho_plane"]==1.0)
        rho_plane(location,holder);
      return;
    }

    void SoundWave::rho_rectangle(double* location, double* holder){
      int dim = (int)para["dim"];
      int nvar = (int)para["nvar"];
      double cen[dim], wid[dim];
      std::cout << "rho_rectangle here" << std::endl;
      if(dim>=1){
        cen[0]=para["cen1"];
        wid[0]=para["wid1"];
      }
      if(dim>=2){
        cen[1]=para["cen2"];
        wid[1]=para["wid2"];
      }
      if(dim==3){
        cen[2]=para["cen3"];
        wid[2]=para["wid3"];
      }
      std::cout << "cen 1: " << cen[0] << " :2 " << cen[1] << " :3 " << cen[2];
      std::cout << "wid 1: " << wid[0] << " :2 " << wid[1] << " :3 " << wid[2];

      holder[0]=1.0;
      for(int i=0; i<dim;i++)
        holder[0]*=(std::abs(location[i]-cen[i])<wid[i]? 1.0 : 0.0);  
      for(int i=1; i<nvar; i++)
        holder[i]=0.0;
      return;
    }

    void SoundWave::rho_circle(double* location, double* holder){
      int dim = (int)para["dim"];
      int nvar = (int)para["nvar"];
      double cen[dim];
      double radius;
      double distance=0;
      // std::cout << "rho_circle" << std::endl;
      if(dim>=1)
        cen[0]=para["cen1"];
      if(dim>=2)
        cen[1]=para["cen2"];
      if(dim==3)
        cen[2]=para["cen3"];
      radius=para["r"];
      for(int i=0; i<dim;i++)
        distance+=std::pow(location[i]-cen[i],2);   
      distance=std::sqrt(distance);
      // 
      holder[0]=(distance<radius? 1.0 : 0.0);
      for(int i=1; i<nvar; i++)
        holder[i]=0.0;
      return;
    }

    void SoundWave::rho_plane(double* location, double* holder){
      int dim = (int)para["dim"];
      int nvar = (int)para["nvar"];
      double cen[dim], k[dim];
      if(dim>=1){
        cen[0]=para["cen1"];
        k[0]=para["k1"];
      }
      if(dim>=2){
        cen[1]=para["cen2"];
        k[1]=para["k2"];
      }
      if(dim==3){
        cen[2]=para["cen3"];
        k[2]=para["k3"];    
      }
      double dotwithnormal=0;
      for(int i=0; i<dim; i++)
        dotwithnormal=dotwithnormal+k[i]*(location[i]-cen[i]);
      //
      // holder[0]=(dotwithnormal>0? 1.0 : 0.1);    
      holder[0]=(0.1+(1.0+tanh(30.0*dotwithnormal))/2.0)/1.1;
      for(int i=1; i<nvar; i++)
        holder[i]=0.0;
      return;
    }
 } // Initial
}  // Dune
